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CLc ABSTRACT 

2 . Context. LS I +61 303 is a puzzling Be/X-ray binary with variable gamma-ray emission up to TeV energies. The nature of the 
' compact object and the origin of the high-energy emission are unclear. One family of models invokes particle acceleration in 
^ shocks from the collision between the B-star wind and a relativistic pulsar wind, whereas another centers on a relativistic jet 
' ' . powered by accretion from the Be star decretion disc onto a black hole. Recent high-resolution radio observations showing a 
putative "cometary tail" pointing away from the Be star near periastron have been cited as support for the pulsar-wind model. 
Aims. We wish to carry out a quantitative assessment of these competing models. 

Methods. We apply a "Smoothed Particle Hydrodynamics" (SPH) code in 3D dynamical simulations for both the pulsar-wind- 
' interaction and accretion-jet models. The former yields a dynamical description of the shape of the wind-wind interaction 
, surface. The latter provides a dynamical estimation of the accretion rate under a variety of conditions, and how this varies with 
orbital phase. 

Results. The results allow critical evaluation of how the two distinct models confront the data in various wavebands. When one 
accounts for the 3D dynamical wind interaction under realistic constraints for the relative strength of the B-star and pulsar 
' winds, the resulting form of the interaction front does not match the putative "cometary tail" claimed from radio observations. 
, On the other hand, dynamical simulations of the accretion-jet model indicate that the orbital phase variation of accretion power 
• • ■ includes a secondary broad peak well away from periastron, thus providing a plausible way to explain the observed TeV gamma 
, ' ray emission toward apastron. 

, Conclusions. Contrary to previous claims, the colliding-wind model is not clearly established for LS I -1-61 303, whereas the 
5_j ■ accretion-jet model can reproduce many key characteristics, such as required energy budget, lightcurve, and spectrum of the 
^ observed TeV gamma-ray emission. 
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1. Introduction 

Some massive X-ray binaries emit electromagnetic radia- 
tion across the entire spectrum, from radio to TeV gamma- 
rays. Cherenkov gamma-ray telescopes have recently de- 
tected variable emission from the binary systems PSR 
B1259-63 (Aharonian et al. 2005a), LS 5039 (Aharonian 
et al. 2005b, 2006a), LS I +61 303 (Albert et al. 2006) and 
Cygnus X-1 (Albert et al. 2007). The first source is a radio 
pulsar in a long and eccentric orbit around a Be star. Its 
high-energy emission is interpreted as inverse Compton 
(IC) upscattering of stellar photons by relativistic elec- 
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trons accelerated in a colliding wind region (e.g. Tavani 
& Arons 1997, Kirk et al 1999, Khangulyan et al. 2007). 
The case of the next two sources is not so clear, since the 
nature of the compact object has not been establishe43 
and jet-like radio features have been detected (Paredes et 
al. 2000, Massi et al. 2001). Nonetheless, some authors 
(e.g. Dubus 2006a) have suggested that all these sources 
may be similar to PSR B1259-63, in the sense of being 
powered by colliding winds. In contrast, Cygnus X-1 is a 
well-established microquasar. 

In the case of LS I +61 303, recent radio imaging 
observations along the orbital period have been cited as 
evidence that the compact object is an energetic pulsar 



^ We note that Casares et al. (2005a) suggest that LS 5039 
might be a black hole. 
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around which a synchrotron "cometary tail" is formed 
(Dhawan et al. 2006). Colliding wind models, however, 
face many problems for explaining the spectral energy 
distribution (SED), the time evolution of the emission, 
the required energy budget, or even the radio morphol- 
ogy. Alternative accretion models, wherein the compact 
object is a black hole and the non-thermal emission arises 
from relativistic jet flows, remain an open possibility (e.g. 
Bosch- Ramon et al. 2006, Orellana & Romero 2007). 

In this paper we apply a Smoothed Particle 
Hydrodynamics (SPH) code to develop more realistic, 
three-dimensional (3D), dynamical simulations of both 
types of model, and use the results of these simulations 
to evaluate the respective merits of each model for match- 
ing observational characteristics of LS I -|-61 303. For the 
pulsar-wind-collision model a key aim is to compute the 
3D form of the wind-wind interaction front under realistic 
constraints for the relative strengths of the B-star and pul- 
sar winds. We adopt a non-relativistic pulsar wind model 
here, but with a wind momentum adjusted to represent a 
relativistic wind with the required total energy. 

For accretion models, a key point to establish is the 
level and phase variation of accretion power in this sys- 
tem. Until now, accretion models for LS I -1-61 303 have 
applied simplistic estimates based on Bondi-Hoyle scaling 
laws (e.g. Marti & Paredes 1995, Gregory & Neish 2002), 
ignoring the perturbation of the compact object on the 
decretion disc of the Be star, as well as the details of the 
accretion disc formation. Our 3D dynamical simulations of 
the accretion disc include cases where the compact object 
is either a neutron star or a black hole, allowing a quanti- 
tative assessment of how well the accretion can be used to 
power jets that might reproduce the observed spectral en- 
ergy distribution and its variation with orbital phase. We 
can then compare the accretion scenario with the collid- 
ing wind one, discussing the critical issues for each model 
type. 

The structure of the paper is as follows. In the next 
section we review the main known properties of LS I 
-|-61 303. Then we outline the main features of the ac- 
cretion/ejection scenario for this source. 

In Sections [4] and [5] we respectively describe the simu- 
lations for the accretion and colliding wind scenarios and 
present the main results derived from them. Section [6] pro- 
vides a global assessment of both competing scenarios and 
a brief summary. 

2. The puzzling gamma-ray binary LS I +61 303 

LS 1+61 303 is a radio-emitting Be X-ray binary discov- 
ered by Gregory & Taylor (1978). The primary star in the 
system is a BO-BO. 5Ve star (Paredes & Figueras 1986), 
with a dense equatorial mass ejection that forms a circum- 
stellar disc (CD). The distance to the system is ~ 2 kpc 
(Frail & Hjellming 1991, Steele et al. 1998). The orbital pe- 
riod is well known from radio observations: P — 26.4960 
days (Gregory 2002). The orbital parameters have been 
determined by Casares et al. (2005b) and Grundstrom et 



al. (2007). The orbit is inferred to be quite eccentric, with 
eccentricity estimates from e — 0.72 ± 0.15 (Casares et 
al. 2005b) to e = 0.55 ± 0.05 (Grundstrom et al. 2007). 
Inferred masses are '--^ 12 Mq for the primary and ~ 2.5 
Mq for the compact object, assuming an inclination angle 
i = 30°, but lower masses for the secondary are possible 
for larger inclinations (e.g. Dubus 2006a). The luminosity 
of the Be star is estimated at ~ 10'^* erg s~^ (Hutchings 
& Crampton 1981). 

Soon after its discovery, LS H-61 303 was associated 
with gamma-ray detections, first with the COS-B source 
2EG 135-1-01 (Gregory & Taylor 1978), and later with the 
EGRET source 3EG J0241-h6103 (Kniffen et al. 1997). In 
2006 it was detected by the MAGIC telescope as a variable 
gamma-ray source at high energies {E > 200 GeV, see 
Albert et al. 2006). The maximum in the gamma-ray light 
curve occurs not during periastron passage (at phase (p = 
0.23), but at phase cf) = 0.5 — 0.6. The radio emission is 
stronger between phases ~ 0.45 — 0.9. 

In X-rays the source has been observed by a number of 
instruments including ROSAT, RXTE, ASCA, Chandra, 
XMM-Newton, and INTEGRAL (e.g. Taylor et al. 1996, 
Leahy et a. 1997, Greiner & Rau 2001, Sidoh et al. 2006, 
Chernyakova et al. 2006, Paredes et al. 2007). The X-ray 
flux is variable, changing from high to low states in a few 
hours. The X-ray luminosity, however, is never very high, 
varying in the range 10'^^ — 10'^^ erg s~^ . Sidoli et al. (2006) 
have compiled the whole SED of the source (see their 
Figure 7). No thermal feature is observed at X-rays. The 
peak of non-thermal emission is in gamma-rays, reaching 
a luminosity of ^ 10^^ erg s~^ in EGRET's energy range. 

At E > 200 GeV the luminosity is - lO^"* erg s'^. All 
this shows that extremely energetic particles are present 
and that a very large energy budget is required to power 
the source, i.e. L > 10'^^ — 10'^'' erg s~^. 

3. The accretion/ejection scenario for 
LS I +61 303 

The idea that the non-thermal emission from LS I -f 61 303 
is powered by accretion onto a compact object was first 
proposed by Taylor & Gregory (1982, 1984). When Massi 
et al. (2001, 2004) detected jet-like radio features extend- 
ing to ^ 400 AU from the core, the source was classified as 
a microquasar (MQ). Since the mass of the compact ob- 
ject is not well-constrained, several authors have proposed 
that the secondary could be a low mass (~ 2.5 Mq) black 
hole (Punsly 1999, Massi 2004, Bosch-Ramon et al. 2006, 
Orellana & Romero 2007). The scenario of an accreting 
pulsar seems to be ruled out by the low-level of the X-ray 
luminositjd, unless most of the accreting matter is ejected 
by the pulsar (Zamanov 1995, Romero et al. 2005). 

If the power engine of LS I -f 61 303 is accretion onto 
a black hole, coUimated outflows of relativistic particles 

^ In the case of an accreting pulsar, the infalling matter 
should impact onto the solid surface, producing strong and 
hard X-ray emission, which is not observed. 
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can be expected, as in other microquasars and in extra- 
galactic quasars. The usual assumption of accreting mod- 
els is that the accretion rate is coupled to the jet power: 
Ljet = Qjct-^accrC^j where gjet < 1 and could be as high as 
0.1. The accretion/ejection models can be divided in two 
groups according to how the high-energy emission orig- 
inates: leptonic models (e.g. Bosch-Ramon et al. 2006, 
Bednarek 2006, Gupta & Bottcher 2006), wherein the 
gamma-rays are the result of IC interactions of relativistic 
electrons in the outflow with photons from the star, the 
CD or locally produced by synchrotron mechanism; and 
hadronic models (e.g. Romero et al. 2003, 2005; Orellana 
& Romero 2007), wherein the emission at £^ > 1 GeV is 
dominated by gamma rays from the decay of neutral pi- 
ous that originate in inelastic pp interactions. Both types 
of models predict variability over the orbital period, since 
one expects orbital phase variation in both the power of 
the jet, as well as in the target density for the relativis- 
tic particles (either photons or ambient matter). An ad- 
ditional source of variability is the absorption of gamma- 
rays in the photon fields of the primary star and the CD 
(e.g. Bosch- Ramon et al. 2006), which can lead to the for- 
mation of electromagnetic cascades and the suppression of 
the gamma-ray signal during periastron passage (Orellana 
& Romero 2007). 

It may be possible in principle to differentiate be- 
tween leptonic and hadronic models through the cutoff 
of the high-energy emission, since electrons cool very fast 
through IC losses in the ambient photon fields, so their 
highest energy can hardly go beyond a few TeV. This is 
not the case for protons, which can reach higher energies. 
In addition, neutrino emission might be detectable from 
this type of sources if the bulk of the gamma-ray emis- 
sion is of hadronic origin (e.g. Christiansen et al. 2006, 
Aharonian et al. 2006b). 

An obvious requirement of the accretion/ejection mod- 
els for the gamma-ray production in LS I -1-61 303 is that 
the accretion rate must be high enough to sustain the 
gamma-ray production with luminosities of ~ 10'^'* erg 
s~-^ at phase </> = 0.5 — 0.6, where MAGIC detected the 
source. We note that the detected EGRET flux, which is 
probably due to LS I -|-61 303, imposes even more severe 
constraints (L^ ^ 10^^ erg s~^). There might be some con- 
tribution from background sources, but most of the fiux 
should come from the binary system, since the emission 
is variable (e.g. Torres et al. 2001, Massi 2004). We will 
discuss this further in Section ID 

4. Numerical simulations of the accretion disc 
formation around the compact object in 
LS I +61 303 

4.1. The code 

Simulations presented here were performed with a three 
dimensional, smoothed particle hydrodynamics (SPH) 
code. The code is basically the same as that used by 
Okazaki et al. (2002) and Hayasaki & Okazaki (2004, 2005, 



2006). It is based on a version originally developed by Benz 
(Benz 1990; Benz et al. 1990) and then by Bate and his 
collaborators (Bate, Bonnell & Price 1995). Using a vari- 
able smoothing length, the SPH equations with a stan- 
dard cubic-spline kernel are integrated with an individual 
time step for each particle. In our code, the Be disc and 
the accretion disc are modeled by an ensemble of gas par- 
ticles with negligible self-gravity, while the Be star and 
the compact object are represented by sink particles with 
the appropriate gravitational mass. Gas particles that fall 
within a specified accretion radius are accreted by the sink 
particle. 

In order to allow the conversion of the kinetic energy 
to heat due to viscosity and to deal with shocks, our code 
uses the standard form of artificial viscosity with two free 
parameters asPH and /3sph, which respectively control the 
strength of the shear and bulk viscosity components and 
that of a second-order, von Neumann-Richtmyer-type vis- 
cosity (Richtmyer & Morton 1994). 

4.2. The simulations 

To optimize the resolution and computational efficiency of 
our simulations, we elected to break the decretion and ac- 
cretion portions of the computation into two separate, but 
linked parts. The first focuses on the decretion of the Be 
disc, starting from the stellar equatorial surface radius i?^, 
then following its evolution under the combined gravita- 
tional infiuence of the B-star and compact object. Because 
of the strong photoionization heating from the stellar ra- 
diation, this Be disc is taken to be isothermal at 0.6r^ 
(Carciofi & Bjorkman 2006), where T^, is the effective tem- 
perature of the Be star. The outward viscous diffusion of 
material in this Be disc provides an effective source for 
capture onto the accretion disc around the compact ob- 
ject, taken to occur whenever material reaches a variable 
capture radius set by 0.9rL, where tl is the Roche-lobe 
radius for a circular binary given approximately by (e.g., 
Paczynski 1971, Warner 1995) 

/ \ 1/3 

^ 0.462 (^^j D. (1) 

Here the mass ratio q = Mx/M^,, where Mx and are 
respectively the masses of the compact object and Be star, 
and D is the instantaneous distance between the stars. 

Following the approach adopted by Hayasaki & 
Okazaki (2004), the second part of the simulation focuses 
on the accretion disc around the compact object, with 
the above phase-dependent mass transfer from the Be 
disc simulation now used as an outer boundary condition 
for the accretion process. The inner edge of the accre- 
tion disc is set to a small fraction of the semi-major axis, 
2.5 X 10~'^a, below which the further accretion onto the 
compact object is expected to occur on rapid time scale 
that is much shorter than that for the simulated region. 
The value of the semi-major axis is a = 6.4 x 10^^ cm. 
Since the role of external heating is less clear than for the 
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Be decretion disc, the material in this accretion disc is as- 
sumed to follow a polytropic relation with index 7 = 1.2, 
representing a compromise between the isothermal (7=1) 
and adiabatic limits (7 = 5/3). 

The numerical viscosity is adjusted to keep the 
Shakura-Sunyaev viscosity parameter ass constant, us- 
ing the approximate relation ass = 0.1 asPH h/H and 
/3sPH = 0, where h and H are respectively the smooth- 
ing length of individual particles and the scale-height of 
the decretion/accrction disc (for additional details, see 
Okazaki et al. 2002). In the Be decretion disc, this scale 
is related to the distance r from the Be star by: 



r \Vait) V-R* 



1/2 



(2) 



where Cg is the isothermal sound speed and Vcmt = 
{GMi,/ Ri.Y/'^ is the critical velocity of the Be star. 

We set the binary orbit in the x-y plane with the ma- 
jor axis along the x-axis. In Be disc simulations, the mass 
ejection mechanism from the Be star is modeled by con- 
stant injection of gas particles at a radius just outside 
the equatorial surface. For the primary, we assume a BOV 
star of = 12Mq, = lOi?©, and = 26, 000 K. 
The scale-height of the CD is ~ 0.03i?^ at r = i?*, then 
increasing outward as r"^/^. For the companion, we take 
Mx = 2.5Mq. The base density of the Be disc is taken 
to be 5 X 10"^^ gcm~'^. Note that the Be star has a two- 
component extended atmosphere: a polar, low-density and 
fast wind 10'^ km s~^) and a cool, slow CD. The matter 
supplied from the equatorial stellar surface drifts outwards 
because of viscous diffusion. The velocity of the CD out- 
flow, within lOi?^, is less than a few km s~^. 

For the binary, we adopt an orbital period Porb = 
26.496 d and an eccentricity e — 0.72 (Casares et al. 
2005b). 
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4.3. Results 

Let us now examine results from a set of simulations with 
the Be disc plane set coplanar with the binary orbital 
plane. Figure [1] shows snapshots around periastron pas- 
sage, which takes place at the orbital phase 4> = 0.23, 
from (a) the Be disc simulation and (b) the correspond- 
ing accretion disc simulation. These snapshots are taken 
after the Be/accretion disc begins to show a regular or- 
bital modulation. Each panel shows the surface density 
over four orders of magnitude, with a logarithmic scaling 
for the greyscale. In Figllji, where the centre of mass of 
the system is fixed, the dark spot near the centre is the 
Be star, whereas the other dark spot represents the com- 
pact object with the variable accretion radius. In Fig[T)3, 
the dark spot at the centre is the compact object with the 
fixed accretion radius of 2.5 x 10~'^a. 

Figd] illustrates the interaction between the Be disc 
and the compact object. In this highly eccentric system, 
strong interaction occurs only for a very short span of 
time around periastron passage. At periastron, the com- 



Fig. 1. Snapshots from (a) a Be disc simulation and (b) 
the corresponding accretion disc simulation and (c) their 
mosaic, covering ^ 0.15 Porb around periastron passage 
at = 0.23. Each panel shows the surface density in cgs 
units on a logarithmic scale. The annotations in each panel 
give the orbital phase and total number of SPH particles. 
Indicated in the first panel are the size and surface density 
scales. 

pact object tidally deforms the Be disc and captures parti- 
cles in the outermost part (Fig[l}i). Matter infalling from 
the Be disc excites the m = 1 density wave in the accre- 
tion disc (FigdiD). At the same time, the accretion disc as a 
whole shrinks by the tidal torque from the Be star. In later 
phases, the density wave propagates inward, enhancing the 
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Fig. 2. Phase dependence of the accretion rate onto the 
compact object. To reduce the fluctuation noise, the data 
is averaged over three orbital period. The vertical dashed 
lines denote periastron passage at phase 0.23. 

accretion rate, as found by Hayasaki & Okazaki (2005). 
The accretion disc gradually expands until the next peri- 
astron passage of the Be star. 

Figure O shows the phase dependence of the accretion 
rate onto the compact object, which is averaged over three 
orbital periods to reduce the fluctuation noise. Note that 
this accretion rate has two peaks per orbit, i.e., a narrow 
peak at periastron, and a broad peak at a later phase. The 
periastron peak (M ^ 5.2 x 10^^ gs~^) is due to shrinking 
of the accretion disc. In our simulations, particles crossing 
the accretion radius of the compact object are taken as 
having accreted and so are removed from the simulation. 
Since the accretion disc in the real system should have a 
much smaller inner radius, this peak could be artificial. To 
test this, we carried out a simulation with the accretion 
radius doubled, to 5 x 10^'^a. We find that the periastron 
peak is indeed higher, about five times the value of the 
standard case, but the secondary, broad peak has a simi- 
lar amplitude in both simulations. As such, it seems that 
periastron peak amplitude may be an artifact of the lim- 
ited resolution, but the extension over a broad secondary 
peak seems more robust. 

This secondary peak has value of M 5.2 x 10^^ gs~^ 
and occurs at a phase ~ 0.52 that lags the periastron pas- 
sage by about 0.3 in phase. As noted above, this secondary 
peak stems from the inward-propagating, m = 1 density 
wave in the accretion disc, which is excited by the ram 
pressure of material transferred from the Be disc at peri- 
astron. The phase lag of ^ 0.3 results from a combination 
of the time for wave excitation and the wave propaga- 
tion time (for a more detailed discussion, see Hayasaki & 
Okazaki 2005). 

These results have important implications for accre- 
tion models of LS I +61 303, since they indicate that a 
peak in non-thermal emission could occur well past pe- 
riastron passage, after phase 0.5, much as observed by 
MAGIC (Albert et al. 2006). Moreover, the energy avail- 
able is more than adequate to explain the observed high- 
energy emission. At phase ~ 0.5 the accretion power is 
~ 4.7 X 10'^'' erg s~^. Assuming that 10% of this power 



goes into an outflow - with just 10% of that converted 
to relativistic protons -, then a gamma-ray luminosity 
of roughly ^ 8 x 10^"* erg s~^ is expected from inelastic 
pp interactions with the medium surrounding the com- 
pact object (Orellana & Romero 2007). If the radiative 
properties of the jet are dominated by relativistic elec- 
trons or electron/positron pairs, gamma-ray production 
through IC interactions could also explain the observed 
emission with the inferred accretion rates (Bosch-Ramon 
et al. 2006). Since our simulations suggest the periastron 
peak at phase 0.23 may be artificial, opacity effects do not 
need to be invoked to explain the near-apastron phase of 
gamma-rays observed by MAGIC. However, as discussed 
in Section [HI in practice such opacity effects may indeed 
be significant, making it difficult to detect TeV gamma 
rays emitted near periastron. 

5. Colliding winds and cometary tails 

The colliding wind scenario for LS I +61 303 was first 
proposed by Maraschi & Treves (1981), and then revis- 
ited by Dubus (2006a). Recently, additional authors have 
proposed that an energetic pulsar with a relativistic wind 
could power the gamma-ray source (Chernyakova et al. 
2006, Neronov & Chernyakova 2007). The central idea 
behind these models is that the collision between the 
relativistic pulsar wind and the slow wind from the Be 
star should result in a stong shock, thus providing a site 
for particle acceleration, much as is thought to occur in 
other massive binaries like WR 140 (Eichler & Usov 1993, 
Benaglia & Romero 2003). Relativistic electrons under- 
going IC interactions with stellar photons then produce 
the observed gamma-ray emission. Since these electrons 
should also produce synchrotron radiation in any ambient 
magnetic fields, the colliding wind region itself could be 
detectable as an extended radio source, perhaps even with 
an elongation directed away from the Be star (Mirabel 
2006, see Figure 1, right panel, and Dubus 2006a). 

Dhawan et al. (2006) present a series of VLBA im- 
ages, spaced 3 days apart, and covering the entire orbit 
of LS I +61 303. The images show extended radio struc- 
tures, with the phase near periastron having indeed an 
apparent elongation away from the primary star. Since no 
macroscopic relativistic velocities are measured, the au- 
thors identify this elongation as "a pulsar wind nebula 
shaped by the anisotropic environment, not a jet". Citing 
the known pulsar PSR B1259-63 as example, they con- 
clude that LS I +61 303 is a Be-pulsar binary, and not a 
M^. 

However, there are several crucial issues for this inter- 
pretation. First, it is important to note that, unlike PSR 
B1259-63, no pulsation is detected in LS I +61 303. As 
such the viability of such a pulsar model requires either 
that any pulsed signals must be quenched by a sufficient 
column of matter throughout the orbital phase (including 

Note, however tfiat no extended radio emission have been 
found so far in PSR B1259-63. 
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at apastron), or, alternatively, that there is a special ge- 
ometry between the magnetic and rotational axis of the 
putative neutron star. 

Second, for a reasonable gamma-ray production efh- 
ciency of a few percent, the gamma-ray luminosity inferred 
by EGRET, lO'^^ erg s~^, requires the total power of 
the pulsar wind to be higher than several times 10'^^ erg 
s~ . Assuming a power of 10^^ erg s-\ Dubus (2006a) 
cannot reproduce the observed SED, as seen through his 
Figure 8, which shows the EGRET flux to be an order of 
magnitude above the predicted emission. 

Finally, this large required pulsar wind power implies 
also a large pulsar wind momentum, requiring then an 
even stronger Be-wind momentum if it is to produce the 
elongated, cometary shape inferred in the radio images. 
The key parameter determining the shape of the interac- 
tion surface of colliding winds is the ratio of wind momen- 
tum fluxes, given by 



where Vbc and Mbc are the velocity and mass loss rate 
of the Be wind, and -EpsR is the power of the pulsar 
wind. Taking Vbc = lO^kms"^ and Mbc = IO-^Mqjt-^ 
as strong wind parameters for a BOV star, we have rj ^ 
0.53^;psR/1036ergs-^ In a simple 2D model that Ignores 
orbital motion, the interaction surface approaches a cone 
with opening half-angle, 

= 180° , (4) 

which can be derived from the analysis of Antokhin et al. 
(2004). Note then that for the above strong Be-wind pa- 
rameters, even the minimally required pulsar wind with 
energy ii^psR ~ 10'^^ erg s^^ would force a quite wide in- 
teraction cone, with half-angle 9 « 62°. And for somewhat 
larger, but still quite modest pulsar wind energies EpsR > 
2 X lO^^ergs"^, the pulsar momentum overwhelms the 
Be-wind, forcing an opening half- angle 9 > 90°, implying 
that the interaction cone is actually wrapped around the 
Be-star, not the pulsar. (For more detailed discussion of 
the shape of the interaction surface, see, e.g., Antokhin 
et al. 2004). Thus even in such simplified 2D interaction 
models, it is problematic to obtain an interaction surface 
that is elongated away from the Be-star, as inferred in the 
radio images. 

The situation is further complicated if one takes into 
account the azimuthal deflection associated with orbital 
motion, which can be particularly significant near perias- 
tron in this highly eccentric system. To explore how such 
orbital effects alter the shape of wind interaction surface, 
we have carried 3D SPH simulations of the collision be- 
tween a Be stellar wind (with the above standard mass loss 
rate and fiow speecQ) and a much faster flow intended to 

* We remind the reader that the relevant Be-wind here is the 
fast stellar wind with spherically symmetric geometry and not 
the equatorial CD. 



represent the pulsar wind. Since our SPH simulations can- 
not directly simulate a relativistic fiow, we have adopted 
a fast but non-relativitic wind, with a speed of 10"* kms~^ 
for the pulsar wind, adjusting the mass-loss rate so as to 
give the same momentum as a relativistic flow with the 
assumed energy, thus yielding a momentum ratio given as 
above by 77 ~ 0.53-EpsR/10^^ergs~^. Our approach is de- 
signed to account for the relative ram pressure balance be- 
tween the colliding winds, since this is the primary factor 
in defining the shape of the wind-wind interface. However, 
by assuming a cooled post-shock flow we are neglecting the 
effect of the gas pressure. We are also neglecting the effect 
of the Be disc, because it is unlikely that a geometrically 
thin equatorial outflow with very low velocity could pro- 
duce the broad "cometary tail" inferred from the radio 
maps. In the future we plan to investigate models with a 
more complete energy balance treatment, including shock 
heating and radiative and other relevant cooling processes. 

For simplicity, we take both winds to be isothermal 
and coasting without any net external forces, assuming 
in effect that gravitational forces are either negligible (i.e. 
for the pulsar wind) or are effectively canceled by radiative 
driving terms (i.e. for the Be star). As in the Be and ac- 
cretion disc simulations, we also ignore the self-gravity of 
the wind flows. The artificial viscosity parameters adopted 
are aspH = 1 and /3sph = 2, since we deal now with a sit- 
uation with no shear and a lot of compression. We set the 
binary orbit in the x-y plane and the major axis of the or- 
bit along the x-axis. The outer simulation boundary is set 
at r = 5.25a from the centre of mass of the system. This 
means scales of ~ 1 mas as seen from the Earth. Particles 
crossing this boundary are removed from the simulation. 

Figure [3] shows a snapshot of a colliding wind model 
for LS I +61 303 with i^PSR = 1036ergs-\ and thus 77 = 
0.53. On a logarithmic scale with cgs units, the greyscale 
shows the density in (a) the orbital plane and (b) the 
plane perpendicular to this orbital plane and through the 
orbit's major axis. Darker and lighter regions represent 
respectively the Be wind and the pulsar wind. The dark 
spot near the origin of the figure represents the Be star, 
and a small dark spot on the left side of the Be star is 
the pulsar. Note that the density in the pulsar wind is 
negligible in this simulation. 

Although the interaction surface in the simulation ex- 
hibits variations from instabilities, its global shape is eas- 
ily traced and illustrates clearly our key central point. 
Namely, when orbital effects are included, even the most 
favorable assumptions toward a large Be/pulsar wind mo- 
mentum ratio do not produce the simple elongated shape 
inferred in the VLBI radio image, which was previously 
cited as strong evidence in favor of a pulsar wind interac- 
tion scenario. 

6. Discussion and conclusions 

The previous section shows that the colliding wind model 
has serious problems, both in matching the observed mor- 
phology at radio wavelengths, and in accounting for the 
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Fig. 3. Wind collision interface geometry at periastron 
for collision between a Be-star wind with a pulsar wind of 
energy Epsr — lO'^^ergs"^, corresponding to pulsar/Be 
wind momentum ratio of = 0.53. The greyscale plot 
shows the density in (a) the orbital plane and (b) the plane 
which is perpendicular to the orbital plane and through 
the major axis of the orbit, on a logarithmic scale with cgs 
units. The dark spot near the origin represents the Be star, 
while the small dark spot to the left represents the pulsar. 
The inset shows the close-up view of the region near the 
stars. Annotations in the figure give the orbital phase 
and the number of particles in the Be wind, A^i , and in the 
pulsar wind, The outer simulation boundary is set at 
a radius r = 5.25a from the centre of mass of the system. 



observed gamma-ray emission through a pulsar wind with 
a power of the order of 10^^ erg s~^. An additional problem 
is the observed gamma-ray phase variation, which shows a 
peak after phase 0.5. Detailed modeling of adiabatic losses 
of electrons or very strong Compton losses during perias- 
tron passage (resulting in a high-energy cutoff) could ex- 
plain a low flux at phase 0.23 (see Khangulyan et al. 2007 
for a treatment of this kind in the case of PSR B1259- 
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63) ; but it remains difficult to explain why the gamma-ray 
source is not detected before periastron passage, as it is 
in PSR B1259-63. 

Orientation effects might be invoked to explain such a 
situation, but detailed modeling is lacking. In the accre- 
tion models, these problems are solved in a natural way 
by the variable accretion rate (see Section d]). 

Opacity effects to gamma-ray propagation in the ambi- 
ent photon fields of the Be star and disc offer an additional 
possibility for quenching the high-energy emission close to 
periastron passage. Our SPH simulations provide the az- 
imuthally averaged surface density profiles for the Be disc, 
and together with the assumption that the disc is isother- 
mal with a temperature Tdisc = 15, 600 K, this allows us 
to model the emission of disc photons. Adding then the 
photon field of the Be star when accounting for the finite 
stellar size (Dubus 2006b), we apply the basic approach of 
Becker & Kafatos (1995) to compute the photon-photon 
attenuation of gamma-rays generated at the base of the 
jet (a height ~ 10 cm above the black hole), and di- 
rected toward the observer (see Appendix A for further 
details). Figure |4] plots the phase variation of the total 
opacity along the orbit, for two different photon energies. 
We see that at periastron the TeV gamma-rays have an 
optical depth r w 2, corresponding to a reduction by a 
factor exp(— 2) « 0.14, implying that longer integrations 
will be required to detect the attenuated source. But since 
the absorption is low during the rest of the orbit, collid- 
ing wind models might predict a detection not only after, 
but also before periastron passage, unless effects related 
to the viewing angle are invoked. In any case, the collid- 
ing wind model would also face some problems to explain 
the lightcurves at other wavelengths, where the radiation 
peaks at similar phases as the VHE emission. 

We note that the accretion model is also not free of 
problems. The changing morphology of the radio emis- 
sion (Dhawan et al. 2006) requires a highly unstable jet. 
This was already suggested by the Massi et al. (2004) ra- 
dio maps. The origin of this instability may be related 
to the density waves in the accretion disc found in this 
paper, but detailed simulations are required to reproduce 
the observed phenomenology. Alternatively, the impact of 
wind inhomogeneities could cause the deflection of the jet 
(Bosch- Ramon et al. 2006). 

Another issue concerning accretion models is the low 
X-ray luminosity observed from this system. The peak 
accretion rate, when expressed in Eddington units, is 
m 0.01. Advection-dominated accretion solutions exist 
for the inner accretion disc up to to < TOdit ^ 0.05 — 0.1 
(Narayan et al. 1998, see p. 160). Since the viscosity param- 
eter and TOcrit are related through merit ~ ct^, the models 
discussed in Section 2] are consistent with advective, un- 
derluminous accretion, which could explain the paucity of 
thermal X-ray emission in LS I -1-61 303. 

In addition, since the global change of the accretion 
rate predicted by the simulations is ~ 30 % and the opac- 
ity is negligible close to the apastron, accretion-jet mod- 
els should introduce effects of the changing environment 
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Fig. 4. Optical depth for the propagation of 7-ray photons 
with energy E = 200 GeV (lowest energy of MAGIC de- 
tection) in the direction to the observer, and E = 1 TeV. 
Circumstellar disc and stellar contributions are shown sep- 
arately. The vertical line indicates the phase of periastron 
passage. 

(photon fields for IC scattering, surrounding matter for pp 
interactions) to explain a non-detection at phases beyond 
~ 0.7 through longer exposures of current Cherenkov tele- 
scopes. In this sense, the study of Be-disc mass not cap- 
tured by the compact object is quite interesting since it 
could be used to explain the infra-red excess observed in 
the source (Marti & Paredes 1995) and provide an ad- 
ditional source of perturbations and interactions for the 
jet. 

Rapid changes in the jet should produce rapid vari- 
ability in the gamma-ray emission, that could be observed 
with GLAST. AGILE and GLAST observations are cru- 
cial to determine the MeV-GeV source luminosity and 
spectrum. MAGIC II could be able to detect the source 
at the periastron with significant integrations, yielding in- 
formation on the opacity effects and their origin. Neutrino 
detection or non-detection with ICECUBE will shed light 
on the nature of the gamma-ray emission. 

We conclude that the current evidence favors an ac- 
cretion/ejection scenario for LS I -1-61 303. The debate, 
however, is not closed and we can expect to learn new 
lessons from this fascinating source in the near future. 
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Appendix A: Calculation of the 7-ray opacity in 
the star+disc photon field 

This appendix briefly summarizes the expressions used for 
the calculation of the optical depth t{E) given in Figure 
[31 We follow the treatment given by Dubus (2006b) and 
Becker & Kafatos (1995). The differential absorption prob- 
ability seen by a 7-ray of energy E located at a position 
P and traveling in the direction e-y through a photon field 
with photons of energy e has two terms: 
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dr-f^ — dedl [ (1 — e,., 
(1 -e-y 



e^)n^(e) aj^dQ^, + 

^^disc)^^disc (e) CT^-ydildisc] (A.l) 

where e.^, edisc and are unit vectors in the photon di- 
rection of propagation; dQ is the solid angle, seen from 
P, of the surface emitting the photons (additional angles, 
not shown in the figure, are defined to parametrize each 
surface, star and disc.) The radiation density, n{e), is in 
photons cm"'^ erg~^ sr~^. The length I of the 7-ray path is 
measured from its emission point S, which is at a distance 
do and forming an angle tpo from the center of the star (see 
a sketch in Figure ETTjl . The interaction cross-section a.y^ 
depends only on the square of the total energy of the two 
interacting photons in the center of mass frame (which, in 
turn, depends on the collision angle between the photons, 
Gould & Schreder, 1967): 

1/2 



eE 



■(1 



>r/disc J 



1 



2/3(/32-2) + (3-/34)ln 



1 + /3' 
1-/3. 



,(A.3) 



where rp — e'^/nicC^ is the classical electron radius. The 
integration over e should be done above the threshold en- 
ergy for pair creation, whereas I is integrated up to infinity 
in order to compute the probability of detection by a dis- 
tant observer. 

We have considered the star and the circumstellar disc 
as black-body emitter^, characterized by T^, and i?*, and 
Tdisc and -Rdisc, respectively. At distances greater than a 
few stellar radii the integration over solid angles reduces 
to the point-source approximation. 

The trajectory considered for the 7-ray photons starts 
at the orbital position of the compact star, assuming that 
the high-energy emission originates close to it. The or- 
bital variation of the opacity through the dependence of 
do and ipo with the orbital phase was calculated as in 
Dubus (2006b). 




Fig. A.l. Scheme of the geometry considered for the 
propagation of a gamma ray emitted at S and absorbed 
at position P. The absorber field is formed by photons of 
the stellar surface and its circumstellar disc. 



(A.2) 



^ We note that the emission of outflowing viscous discs of Be 
stars could be properly modeled as in Porter (1999). However, 
in the present case, in order to account for the luminosity at- 
tributed by Casares et al. (2005b) to the circumstellar disc 
contribution in LSI -1-61 303, we have considered a black-body 
approximation. 



